Excitonic condensation in a symmetric electron-hole bilayer 
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Using Diffusion Monte Carlo simulations we have investigated the ground state of a symmetric 
electron-hole bilayer and determined its phase diagram at T = 0. We find clear evidence of an 
excitonic condensate, whose stability however is affected by in-layer electronic correlation. This 
stabilizes the electron-hole plasma at large values of the density or inter-layer distance, and the 
Wigner crystal at low density and large distance. We have also estimated pair correlation functions 
and low order density matrices, to give a microscopic characterization of correlations, as well as to 
try and estimate the condensate fraction. 

PACS numbers: 71.35.-y, 71.10.-w, 73.21.-b 
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Electron-hole systems have attracted a lot of interest 
over the years Jlj. The Coulomb attraction existing be- 
tween the two kinds of Fcrmions naturally brings about 
pairing, and hence the possibility of a coherent state [Q. 
It was soon realized (| that systems of spatially sepa- 
rated electrons and holes, such as a bilayer, have a num- 
ber of advantages with respect to conventional bulk sam- 
ples |l| in which electrons and holes occupy the same 
region. Thus, while in a homogeneous semiconductor 
the excitonic condensate would be an insulator in 
a bilayer superconductivity is in principle possible ||. 
The interest in such systems has greatly increased in re- 
cent years, due to the increasing ability to manufacture 
high quality semiconductor quantum well (QW) struc- 
tures, where electrons and holes are indeed confined in 
different regions, between which tunneling can be made 
ncglig ible §. Also, in a number of such systems ex- 
perimental evidence of an excitonic condensate has been 
claimed [|[ |[ 0|. However, even for the simplest two 
dimensional (2D) model, i.e., a symmetric bilayer, the 
theory has been concerned so far with mean-field treat- 
ments [|| ||, |§] based on a BCS like wavefunction, with 
one exception [^o| . 

In this Letter we report on the first extensive computer 
simulations of a symmetric electron- hole bilayer (SEHB). 
We have undertaken this study to assess the effect of elec- 
tron correlation on the phase diagram of the SEHB and 
on the existence of a condensate, which we directly char- 
acterize. To perform simulations, we have resorted to 
fixed-node diffusion Monte Carlo (FN-DMC) @ |l§, a 
method which is stable and variational and is known to 
yield extremely accurate results for homogeneous elec- 
tron systems jl3|, [lij. We should stress that our goal 
is to determine the properties of the simplest 2D model, 
i.e., the SEHB, with unprecedented accuracy, also to pro- 
vide a benchmark against which approximate many-body 
treatments may be tested. Therefore, we are not consid- 
ering here effects such as finite layer thickness, inter-layer 
tunneling, or band anisotropy, which may play an im- 



portant role [L5|, [16| in describing more realistic QW 
structures. Also, to limit the computational load, we 
have restricted our study to three phases: the excitonic 
phase (EP), the spin unpolarized two component plasma 
(2CP), and the triangular Wigner crystal (WC). 

In the absence of magnetic fields, the Hamiltonian of 
such an ideal SEHB reads 
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with e the background dielectric constant, d the inter- 
layer distance, and m e = = m* the common effec- 
tive mass of electrons and holes. In the following we use 
Ry* = e 2 /2ea* B as unit energy, a* B — H 2 e/m*e 2 , and r s a B 
as unit length. As is well known the parameter r s , defined 
in terms of the in-layer areal density n by irr 2 a* B 2 = 1/n, 
measures the in-layer coupling strength. In the SEHB 
one is lead to define an additional parameter, measuring 
the importance of the inter-layer coupling, as the ratio 
of the typical inter-layer and in-layer Coulomb energies, 
namely 7 = 1/d. At T = 0, which is the case considered 
here, the model is completely specified by r s and d or 7. 
In FN-DMC [|1J one propagates a trial wavefunction 
in imaginary time, to project out the higher energy 
components and sample the lowest energy state ^0 with 
the nodal structure of ^t- This establishes a correspon- 
dence between nodal structure, i.e., trial function, and 
phase. To study the EP we resort to a BCS-like trial 
function, which is known to provide a good mean-field 
description M ^, |l7j both at high and at low density. 
In practice we set ^>§ p = D^DnJ, where D aa is a 
determinant of pair orbitals <£>(rf CT — r^ a ) and J is a Jas- 
trow factor, accounting for two-body correlations. We 
choose ip(r) as the exact numerical solution of the mean- 
field problem 0, as in selected cases we found that this 
yields a lower variational energy than other choices [G~8|. 
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TABLE I: Energy per particle of various phases of the SEHB, 
in Ry* , according to FN-DMC. All results are extrapolated 
to the thermodynamic limit [B3J . 



r a 


d 


TE 


2CP 


E 


WC 


1 


0.0 


1.69 


-0.833(8) 


-0.808(9) 




2 


0.1 


0.84 


-0.6947(5) 


-0.6976(7) 




2 


0.2 


1.00 


-0.6116(7) 


-0.6006(7) 




2 


0.5 


1.34 


-0.5405(5) 


-0.5260(4) 




5 


0.2 


0.57 


-0.3732(6) 


-0.3822(2) 




5 


0.5 


0.93 


-0.3125(2) 


-0.3104(1) 




5 


1.0 


1.40 


-0.3009(2) 


-0.2987(2) 




10 


0.5 


0.68 


-0.1801(1) 


-0.18172(9) 




10 


1.0 


1.23 


-0.17153(8) 


-0.17085(4) 




10 


1.5 




-0.17035(1) 




-0.16947(1) 


20 0.05 


0.14 


-0.3275(2) 


-0.3288(1) 




20 


0.5 


0.50 


-0.10051(5) 


-0.10227(3) 


-0.10157(2) 


20 


1.0 


1.08 


-0.09304(1) 


-0.09316(2) 


-0.093176(8) 


20 


1.3 


1.38 


-0.09264(2) 


-0.09261(2) 




20 


1.5 


1.51 


-0.09260(2) 


-0.09248(2) 


-0.092533(5) 


20 


3.0 




-0.09245(2) 




-0.092354(6) 


22 


1.0 




-0.08533(2) 




-0.085484(7) 


22 


2.0 




-0.08482(1) 




-0.084770(6) 


22 


3.0 




-0.08478(1) 




-0.084767(5) 


30 


0.5 


0.45 




-0.07192(2) 


-0.07161(2) 


30 


1.0 


1.00 




-0.064403(8) 


-0.064483(3) 


30 


1.5 


1.48 




-0.063833(6) 


-0.063921(3) 


30 


3.0 




-0.06377(7) 




-0.063827(3) 



For the normal phases the trial function is taken as 
* T = D\D^D\D^J, with D a a a Slater determinant of 
one-particle orbitals (plane waves for the 2CP and gaus- 
sians localized at the crystal sitesjl9| for the WC). The 
Jastrow factor J = exp[-(l/2) Y,^y Y! ifl , iu «M,"( r wJ]i 
with the Greek index denoting particle type and spin 
projection, is built using RPA pseudopotentials [^o], ^lf . 

We have performed DMC simulations for systems with 
N = 58 and N — 56 particles per layer, respectively 
for the 2CP and the EP and for the WC, using peri- 
odic boundary conditions. Finite size effects have been 
mitigated as usual [^2] by performing the Ewald sum- 
mation on the infinite periodic replicas of the simulation 
box and using, for the 2CP, numbers of particles corre- 
sponding to closed shells of orbitals. We have also carried 
out variational Monte Carlo (VMC) simulations, for sev- 
eral values of N (ranging up to 114 or 120 depending 
on the phase), to determine the size dependence j23| of 
the energy. Assuming that this is the same for VMC and 
FN-DMC fHH, H , we have obtained the FN-DMC en- 
ergies in the thermodynamic limit, which we report p(| 
in Table || for the three phases that we have studied. 

Use of the energies of Table | yields the phase dia- 
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FIG. 1: Phase diagram of the SEHB according to the FN- 
DMC. Circles, squares and triangles respectively indicate the 
stability of the 2CP, EP,and WC. Dashed, gray and black lines 
give approximate boundaries between pair of phases obtained 
from the crossing of the energies at given r s with varying d. 
The dotted curve reports an estimate of the 2CP-EP bound- 
ary from an approximate theory (BOl. 



gram shown in Fig. |l|. It is evident that correlation has 
here qualitative effects. Whereas the mean-field predicts 
stability of the EP with respect to the 2CP everywhere, 
at high enough density and/or at large enough distance 
FN-DMC predicts the stability of either the 2CP or the 
WC. Naively, one would expect that when the inter-layer 
coupling 7 — 1/d <C 1 each layer should behave as an iso- 
lated 2D electron gas (2DEG) [f27| . This is the case with 
the symmetric electron bilaye r jpo[ p8[ [29[| , for which the 
phase diagram of the 2DEG [|28f| is quickly recovered as 
d exceeds 1, and it appears to be the case also for the 
SEHB, when one keeps in mind that only the unpolar- 
ized 2CP is considered here. The phase diagram turns 
out to be fairly robust. Thus, use of the finite-iV FN- 
DMC energies leaves it essentially unchanged, while use 
of the VMC energies only brings about a minor change 
in the boundary between the WC and the 2CP, which 
moves at about r s = 20 for d > 1.5. 

In Table |] we have also reported the excitonic radius, 
defined by r 2 E = (ip\r 2 \(p)/((p\ip). It is evident, from an 
inspection of the Table, that a correlation exists between 
the stability of the EP and te being smaller than 1, i.e., 
smaller than the characteristic in-layer length scale. This 
points to the importance of in-layer correlations, which 
are neglected in mean-field (|, Also, it turns out that 
in our simulations the exciton is always much smaller 
than the side of the simulation box, \J Nit ~ 13. Thus, 
our description of the EP should not be affected by the 
finite spacing of the energy levels jl0) . 

A peculiar property of superconductors is the off- 
diagonal long-range order (ODLRO) exhibited by the re- 
duced density matrices in the coordinate space represen- 
tation plf. In a Fermion system ODLRO shows up in 
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d=0.2Extr. DMC 
d=0.5 Extr. DMC 
■ - d=0.2 BCS mean-field 
- d=0.5 BCS mean-field 
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FIG. 2: Projected two-body density matrix of the SEHB at 
r s = 5, according to FN-DMC and BCS mean-field, for N = 
58. The full curves are fits to the tails of the simulation data 
(see text and Table |n| for details) . 



TABLE II: Condensate fraction a of the SEHB according 
to FN-DMC (extrapolated estimate) and VMC, for a system 
with N=58. Here d = d-r s is the interlayer distance in units of 
a* B . Also reported are the BCS mean-field prediction and, for 
FN-DMC, the reduced x 2 of the fit yielding the condensate 
fraction (see text). 



r s 


d 


q(BCS) 


a(VMC) 


a(FN-DMC) 


x 2 


2 


0.2 


0.55 


0.187(4) 


0.284(9) 


0.59 


5 


1.0 


0.74 


0.151(2) 


0.215(4) 


0.22 


5 


2.5 


0.48 


0.095(3) 


0.108(5) 


0.34 


20 


1.0 


0.98 


0.027(1) 


0.020(2) 


0.47 



the two-body density matrix 63], with the appearance 
of an eigenvalue which scales with the number of par- 
ticles. In a translational invariant system, such as the 
SEHB in the excitonic phase, ODLRO implies for the 
two-body density matrix, to leading order in N, the fol- 
lowing asymptotic behavior: 



P2(Xe,x^;x e ,x h ) 
|x e - Xfcl, 



aA7*(|x' e -x',J)/(|x e 

*4- x /J~£> |x e -x 



x h |), (2) 



where a < 1 is the condensate fraction and £ is the range 
of the normalized pair amplitude /(|xi — Xa|). In order to 
estimate the condensate fraction it is convenient to resort 
to the projected two-body density matrix (P2BDM) 



— J dx e dx h p 2 (x e +x,x/j +x;x e ,x /l ) 



(3) 



which tends to a in the large x limit, as is immediately 
found combining Eqs. and (||). 

A simple estimator of h(x) is given by 



(4) 




FIG. 3: Spin resolved electron-hole pair correlation function 
of the SEHB at r s = 20, according to FN-DMC. 



with M c the number of particle configurations used and 
R' obtained from the configuration R by rigidly translat- 
ing an electron-hole pair |}2j by x. In practice, for each R 
we generate a few translations x uniformly distributed in 
the simulation box and, for better statistics, we also aver- 
age Eq. (^) over all electron-hole pairs. As Eq. (Q) only 
yields a mixed estimate in DMC, we have also performed 
VMC calculations to get extrapolated estimates Q| ac- 
cording to ft-Extr(x) = 2fe fl Mc(x) - hvMcfa)- 

An illustration of our results for h(x) is given in Fig. 
||. Indeed, h(x) appears to saturate at large x. In the 
absence of exact information on its asymptotic form, we 
have fitted its tail to q + 4/x 2 , x > 5 and estimated the 
condensate fractions reported in Table |l| for the cases 
studied. It is evident that in-layer correlation, which 
is absent in mean field, causes a substantial reduction 
of the condensate fraction — a reduction which becomes 
more pronounced with increasing the in-layer coupling 
r s , at given interlayer distance, and to less extent with 
increasing the interlayer distance, at given r s . This ef- 
fect is particularly strong at large coupling, yielding the 
essential suppression of the condensate at r s — 20 and 
d = 0.05, with a ~ 0.02 whereas in mean field a ~ 1. 

Additional insight into the nature of the phases of 
the SEHB is provided by the pair correlation functions, 
whose features we briefly summarize here. As the two 
layers are brought together from infinity inter-layer in- 
teraction produces a substantial buildup in inter-layer 
correlations and the screening and weakening of intra- 
layer correlations, as for the electron bilayer |jo| p9[ . 
However, while g e h{f) and g e l(r) monotonously increase 
with decreasing d, near the origin, glHr), first develops a 
correlation hole and then by further diminishing d devel- 
ops a peak. This behavior, which becomes particularly 
marked at large coupling as is evident from Fig. |^, might 
be interpreted as a tendency toward biexciton formation 
for large r s and small d. In fact, gH(r) gives correlations 
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between unpaired (f J.) electrons and holes and, indirectly, 
between f and J. excitons. 

To summarize, by performing extensive quantum sim- 
ulations we have shown that correlation has important 
qualitative effects in determining the phase diagram and 
the excitonic condensate of the SEHB, with respect to 
mean-field. Here, we have chosen not to consider spin 
polarized phases nor inhomogeneities of spin or charge, 
such as in density waves and in liquid- vapor coexistence. 
We shall explore some of these interesting phenomena in 
future investigations. 

We acknowledge useful discussions with many col- 
leagues and among the others with D. Neilson, C. Castel- 
lani, R.K. Moudgil. SDP and GS respectively acknowl- 
edge the holding of INFM fellowships while carrying out 
this work and support by MURST through COFIN99. 



* Present address: INFM -Istituto Nazionale di Fisica della 
Materia - Unita di Roma 

1 Present address: Product and Business Development 
I — | Group, Banca IMI, Milano, Italy 

* senatore@ts.infn.it 

[1] T.M.Rice, Solid State Phys. 32, 1 (1977); J.C.Hensel, 

T.G Phillips, and G.A. Thomas, ibid. 32, 88 (1977). 
[2] L.V. Keldish and Y.V. Kopaev, Fiz. Tverd. Tela 6 

2791 (1964) [Sov. Phys. Solid State 6, 2219 (1965)]; 

A.N.Kozlov and L.A. Maximov, Zh. Eksp. Teor. Fiz. 

48, 1184 (1965)[Sov. Phys. JEPT 21, 790 (1965)]; L.V. 

Keldish and A.N.Kozlov, Zh. Eksp. Teor. Fiz. 54, 1978 

(1968)[Sov. Phys. JEPT 27, 521 (1968)]. 
[3] Yu. E. Lozovik and V. I. Yudson, Pis'ma Zh. Eksp. Teor. 

Fiz. 22, 556 (1975)[JETP Lett. 22 , 274 (1975)]; Solid 

State Commun. 19, 391 (1976); Zh. Eksp. Teor. Fiz. 71, 

738 (1976) [Sov. Phys. JETP 44, 389 (1976)]. 
[4] U.Sivan, P.M. Solomon and H. Shtrikman, Phys. Rev. 

Lett 68, 1196 (1992); B. Kane et al., Appl. Phys. Lett. 

65, 3266 (1994). 
[5] L.V. Butov et al. Phys. Rev. Lett. 73, 304 (1994). 
[6] V.B. Timofeev et al. Phys. Rev. B 61, 8420 (2000). 
[7] J. P Cheng et al, Phys. Rev. Lett. 74, 450 (1995) J. 

Kono et al., Phys. Rev. B 55, 1617 (1997); T.P. Marlow 

et al, Phys. Rev. Lett. 82, 2362 (1999). 
[8] Y.E. Lozovik and O.L. Berman, Physica Scrip. 55, 491 

(1997). 

[9] Xuejun Zhu, P.B. Littlewood, S. Hybersten and T. M. 
Rice, Phys. Rev. Lett 74, 1633 (1995),P.B. Littlewood 
and Xuejun Zhu, Phys. Scripta T 68, 56 (1996). 

[10] Antony Chan, PhD Thesis, Indiana University (1996). 

[11] P. J. Reynolds, D. M. Ceperley, B. J. Alder and W.A. 
Lester, J. Chem. Phys. 77, 5593 (1982); see also C. 



J. Umrigar, M.P. Nightngale and K.J Runge, J. Chem. 
Phys. 99, 2865 (1993). 

[12] For a recent survey of QMC methods see, e.g., Quantum 
Monte Carlo Methods in Physics and Chemistry, ed. M.P. 
Nightingale and C.J. Umrigar (Kluwer, Dordrecht, 1999). 

[13] Y. Kwon, D.M Ceperley, and R.M. Martin, Phys. Rev. 
B 53, 7376 (1996); ibid. 58, 6800 (1998). 

[14] D.Varsano, S.Moroni, and G.Senatore, Europhys. Lett. 
53, 348 (2001). 

[15] S.I. Shevchenko, Phys. Rev. Lett. 72, 3242 (1994). 

[16] S.Conti, G. Vignale and A.H. MacDonald, Phys. Rev. B 
57, 6846 (1998). 

[17] C. Comte and P.Nozieres, J. Phys. 43, 1069 (1982); 
P.Nozieres and C. Comte, ibid. 43, 1083 (1982). 

[18] We tried a variational gaussian and the form of Ref. [uj. 

[19] The electron and hole triangular lattices are chosen in 
the energetically favored AA stacking. 

[20] F. Rapisarda, PhD Thesis, Universita di Trieste (1995). 

[21] For the EP J is constructed using the 2CP pseudopo- 
tentials. Depending on coupling the pseudopotentials be- 
tween paired particles or even all interlayer pseudopoten- 
tials may be dropped, to optimize the variational energy. 

[22] D. Ceperley, Phys. Rev. B 18 3126 (1978). 

[23] The functional forms of Ref. E4, which have proven satis- 
factory also for the electron bilayer |2(J], have been used 
for the normal phases and with minor adaptions (like 
E N = E x + b 2 (r s )/N a (r s ) ) also for the EP. 

[24] B. Tanatar and D.M. Ceperley, Phys. Rev. B 39, 5005 
(1989). 

[25] G. Senatore, S. Moroni, and D. Varsano, Solid. St. 

Comm. 119 333 (2001). 
[26] Our energies do not include the phase independent 

Hartree term 2d/r s , accounting for the electric field 

present if the layers are not separately neutralized. 
[27] We only consider 2DEG studies using Slater- Jastrow trial 

functions — the same chosen here for the normal phases. 
[28] F. Rapisarda and G. Senatore, Aust. J. Phys. 49, 161 

(1996). 

[29] G. Senatore, F. Rapisarda and S. Conti, Int. Jour. Mod. 
Phys. B 13 ,479 (1999). 

[30] Lerwen Liu, L. Swierkowsky and D. Neilson, Physica B 
249, 594 (1998). 

[31] C.N. Yang, Rev. Mod. Phys. 34, 694 (1962). 

[32] To reduce size effects only the pair in the simulation box 
is displaced, while its periodic copies are held fixed. We 
have indeed found that this procedure, which is an ob- 
vious extension of the one introduced in Ref. ^ for the 
one-body density matrix, substantially reduces size ef- 
fects in VMC simulations with N = 42, 58, 114. 

[33] WR Magro and D.M. Ceperley Phys. Rev. Lett. 73, 826 
(1994). 

[34] D.M. Ceperley and M.H.Kalos, in Monte Carlo Meth- 
ods in Statistical Physics, K. Binder, ed. (Springer, N.Y., 
1979). 



